Chapter 8 Differential abundance analysis

8.1 CC acclimation vs post2

8.1.1 Enrichment analysis: Ancombc2

8.1.1.1 Structural zeros

# A tibble: 181 × 14
# Rowwise: 
   genome              present average domain phylum class order family genus species completeness contamination length    gc
   <chr>               <chr>     <dbl> <chr>  <chr>  <chr> <chr> <chr>  <chr> <chr>          <dbl>         <dbl>  <dbl> <dbl>
 1 AH1_2nd_11:bin_000… accli     6.79  d__Ba… p__Ve… c__V… o__V… f__Ak… g__A… s__            100            0.12 3.05e6  0.58
 2 AH1_2nd_1:bin_0000… accli     5.96  d__Ba… p__Ba… c__B… o__B… f__Ma… g__O… s__             99.5          0.1  3.92e6  0.44
 3 AH1_2nd_17:bin_000… accli     2.61  d__Ba… p__Ba… c__B… o__B… f__Ri… g__M… s__             98.1          0.06 2.37e6  0.44
 4 AH1_2nd_17:bin_000… accli     2.42  d__Ba… p__Ba… c__B… o__B… f__Ba… g__B… s__             88.0          3.29 4.26e6  0.4 
 5 AH1_2nd_14:bin_000… accli     2.19  d__Ba… p__Ba… c__C… o__L… f__La… g__   s__             97.1          0.82 3.93e6  0.3 
 6 AH1_2nd_5:bin_0000… accli     1.63  d__Ba… p__Ba… c__B… o__B… f__Ma… g__O… s__             99.6          0.04 3.60e6  0.45
 7 AH1_2nd_12:bin_000… accli     1.52  d__Ba… p__Ba… c__B… o__E… f__Co… g__T… s__            100            7.76 3.79e6  0.31
 8 AH1_2nd_8:bin_0000… accli     1.25  d__Ba… p__Ps… c__A… o__R… f__RU… g__   s__             99.6          0    8.93e5  0.36
 9 AH1_2nd_1:bin_0000… accli     1.15  d__Ba… p__Ps… c__A… o__R… f__UB… g__C… s__             98.0          0    9.35e5  0.28
10 AH1_2nd_12:bin_000… accli     0.911 d__Ba… p__Ba… c__C… o__L… f__La… g__C… s__             72.2          2.05 3.51e6  0.42
# ℹ 171 more rows

8.1.1.2 MAG level

8.2 WC acclimation vs post2

8.2.1 Enrichment analysis: Ancombc2

8.2.1.1 Structural zeros

# A tibble: 134 × 14
# Rowwise: 
   genome              present average domain phylum class order family genus species completeness contamination length    gc
   <chr>               <chr>     <dbl> <chr>  <chr>  <chr> <chr> <chr>  <chr> <chr>          <dbl>         <dbl>  <dbl> <dbl>
 1 AH1_2nd_13:bin_000… accli      2.13 d__Ba… p__Ba… c__C… o__L… f__La… g__   s__             97.0          0.58 3.76e6  0.42
 2 LI1_2nd_2:bin_0000… accli      2.07 d__Ba… p__Ba… c__B… o__E… f__Co… g__C… s__             86.9          4.89 3.16e6  0.31
 3 AH1_2nd_12:bin_000… accli      1.99 d__Ba… p__Ba… c__B… o__L… f__St… g__L… s__Lac…         99.9          0.2  2.05e6  0.38
 4 LI1_2nd_3:bin_0000… accli      1.98 d__Ba… p__Ba… c__B… o__S… f__St… g__S… s__Sta…         91.7          3.54 2.96e6  0.32
 5 LI1_2nd_8:bin_0000… accli      1.47 d__Ba… p__Ba… c__C… o__O… f__Ru… g__A… s__             80.1          3.35 2.72e6  0.62
 6 LI1_2nd_2:bin_0000… accli      1.14 d__Ba… p__Ba… c__B… o__R… f__UB… g__C… s__             91.6          1.86 1.13e6  0.28
 7 LI1_2nd_2:bin_0000… accli      1.10 d__Ba… p__Ba… c__C… o__O… f__Ru… g__R… s__             94.7          2.65 2.84e6  0.53
 8 AH1_2nd_6:bin_0000… accli      1.09 d__Ba… p__Ba… c__C… o__L… f__La… g__J… s__             99.8          0.57 3.69e6  0.45
 9 LI1_2nd_1:bin_0000… accli      1.07 d__Ba… p__Ba… c__C… o__T… f__CA… g__R… s__             83.8          7.54 1.93e6  0.3 
10 AH1_2nd_16:bin_000… accli      1.07 d__Ba… p__Ba… c__B… o__B… f__Ri… g__A… s__             93.6          6.68 3.45e6  0.59
# ℹ 124 more rows

8.2.1.2 MAG level

8.3 CI acclimation vs post2

8.3.1 Enrichment analysis: Ancombc2

8.3.1.1 Structural zeros

# A tibble: 195 × 14
# Rowwise: 
   genome              present average domain phylum class order family genus species completeness contamination length    gc
   <chr>               <chr>     <dbl> <chr>  <chr>  <chr> <chr> <chr>  <chr> <chr>          <dbl>         <dbl>  <dbl> <dbl>
 1 AH1_2nd_1:bin_0000… accli     39.3  d__Ba… p__Ba… c__C… o__L… f__La… g__K… s__             98.5          4.38 4.26e6  0.44
 2 LI1_2nd_7:bin_0000… accli      7.37 d__Ba… p__Ba… c__B… o__E… f__Er… g__B… s__             97.9          2.36 2.62e6  0.3 
 3 AH1_2nd_18:bin_000… accli      4.35 d__Ba… p__Ba… c__B… o__B… f__Ta… g__P… s__             98.7          0.01 4.01e6  0.43
 4 AH1_2nd_20:bin_000… accli      2.87 d__Ba… p__Ba… c__C… o__L… f__La… g__R… s__            100.           0.43 3.73e6  0.42
 5 AH1_2nd_5:bin_0000… accli      2.45 d__Ba… p__Ba… c__B… o__B… f__Ba… g__B… s__             86.4          3.24 5.34e6  0.43
 6 AH1_2nd_1:bin_0000… accli      1.71 d__Ba… p__Ba… c__B… o__B… f__Ta… g__P… s__             99.6          1.94 4.53e6  0.51
 7 AH1_2nd_20:bin_000… accli      1.69 d__Ba… p__Ba… c__B… o__L… f__En… g__E… s__             92.6          6.02 4.15e6  0.4 
 8 AH1_2nd_16:bin_000… accli      1.59 d__Ba… p__Ba… c__C… o__L… f__    g__   s__             98.4          0.14 3.87e6  0.28
 9 LI1_2nd_3:bin_0000… accli      1.48 d__Ba… p__Ba… c__B… o__E… f__Er… g__C… s__             99.8          1.24 3.67e6  0.36
10 AH1_2nd_10:bin_000… accli      1.20 d__Ba… p__Ba… c__C… o__L… f__La… g__E… s__             88.9          0.53 3.62e6  0.48
# ℹ 185 more rows

8.3.1.2 MAG level

8.4 post2 CC and CI

8.4.1 Enrichment analysis: Ancombc2

8.4.1.1 Structural zeros

# A tibble: 174 × 14
# Rowwise: 
   genome              present average domain phylum class order family genus species completeness contamination length    gc
   <chr>               <chr>     <dbl> <chr>  <chr>  <chr> <chr> <chr>  <chr> <chr>          <dbl>         <dbl>  <dbl> <dbl>
 1 AH1_2nd_18:bin_000… cc_pos…   36.2  d__Ba… p__Ba… c__B… o__B… f__Ta… g__P… s__             98.7          0.01 4.01e6  0.43
 2 AH1_2nd_15:bin_000… cc_pos…    7.12 d__Ba… p__Cy… c__V… o__G… f__Ga… g__L… s__             99.8          0.87 2.11e6  0.37
 3 LI1_2nd_7:bin_0000… cc_pos…    4.42 d__Ba… p__Ba… c__C… o__C… f__Cl… g__S… s__             98.3          1.37 3.94e6  0.28
 4 AH1_2nd_9:bin_0000… cc_pos…    4.33 d__Ba… p__Ps… c__G… o__B… f__Bu… g__A… s__             72.8          6.01 3.34e6  0.57
 5 AH1_2nd_7:bin_0000… cc_pos…    3.75 d__Ba… p__Ps… c__A… o__R… f__CA… g__C… s__            100            0.49 1.79e6  0.43
 6 AH1_2nd_1:bin_0000… cc_pos…    2.77 d__Ba… p__Ba… c__B… o__B… f__Ta… g__P… s__             99.6          1.94 4.53e6  0.51
 7 AH1_2nd_7:bin_0000… cc_pos…    1.87 d__Ba… p__Ba… c__C… o__L… f__La… g__J… s__            100            2.39 4.97e6  0.43
 8 AH1_2nd_7:bin_0000… cc_pos…    1.49 d__Ba… p__Ba… c__B… o__B… f__Ta… g__P… s__             93.7          2.1  5.26e6  0.43
 9 AH1_2nd_7:bin_0000… cc_pos…    1.47 d__Ba… p__Ba… c__C… o__O… f__Ru… g__   s__             91.7          4.31 2.93e6  0.49
10 AH1_2nd_19:bin_000… cc_pos…    1.35 d__Ba… p__Ba… c__C… o__L… f__La… g__C… s__             92.8          4.49 5.61e6  0.47
# ℹ 164 more rows

8.4.1.2 MAG level

8.5 post2 WC and CI

8.5.1 Enrichment analysis: Ancombc2

####Structural zeros

# A tibble: 161 × 14
# Rowwise: 
   genome              present average domain phylum class order family genus species completeness contamination length    gc
   <chr>               <chr>     <dbl> <chr>  <chr>  <chr> <chr> <chr>  <chr> <chr>          <dbl>         <dbl>  <dbl> <dbl>
 1 AH1_2nd_14:bin_000… ci_pos…   17.8  d__Ba… p__Ba… c__C… o__C… f__UB… g__   s__             94.8          0    2.02e6  0.36
 2 AH1_2nd_1:bin_0000… ci_pos…    4.77 d__Ba… p__Ps… c__A… o__R… f__UB… g__C… s__             98.0          0    9.35e5  0.28
 3 AH1_2nd_12:bin_000… ci_pos…    4.76 d__Ba… p__Ba… c__B… o__M… f__My… g__M… s__             92.6          0.18 1.01e6  0.27
 4 AH1_2nd_8:bin_0000… ci_pos…    4.34 d__Ba… p__Ps… c__A… o__R… f__RU… g__   s__             99.6          0    8.93e5  0.36
 5 AH1_2nd_1:bin_0000… ci_pos…    4.10 d__Ba… p__Ba… c__C… o__L… f__La… g__M… s__             96.6          0.56 3.88e6  0.33
 6 AH1_2nd_19:bin_000… ci_pos…    3.95 d__Ba… p__Ve… c__V… o__O… f__LL… g__   s__             92.8          0    8.89e5  0.53
 7 AH1_2nd_2:bin_0000… ci_pos…    3.09 d__Ba… p__Ve… c__V… o__O… f__LL… g__   s__             93.8          0    8.62e5  0.56
 8 AH1_2nd_6:bin_0000… ci_pos…    2.21 d__Ba… p__Ba… c__B… o__F… f__We… g__M… s__             98.7          0.1  3.38e6  0.36
 9 AH1_2nd_10:bin_000… ci_pos…    1.72 d__Ba… p__De… c__D… o__D… f__De… g__B… s__             97.9          0.61 3.77e6  0.51
10 AH1_2nd_5:bin_0000… ci_pos…    1.67 d__Ba… p__Ba… c__C… o__O… f__Os… g__P… s__             93.4          7.74 3.92e6  0.58
# ℹ 151 more rows

8.5.1.1 MAG level

8.6 Community level plots

accli_post2<-sample_metadata%>% 
  filter(time_point == "1_Acclimation"|time_point == "6_Post-FMT2")
accli_post2$newID<-paste(accli_post2$type, "_", accli_post2$time_point)

GIFTs_functions_community %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(accli_post2, by = join_by(sample == Tube_code)) %>%
  filter(time_point == "1_Acclimation"|time_point == "6_Post-FMT2") %>%
  select(c(1:29, 36,39)) %>%
  pivot_longer(-c(sample,type,time_point),names_to = "trait", values_to = "value") %>%
  mutate(trait = case_when(
    trait %in% GIFT_db3$Code_function ~ GIFT_db3$Function[match(trait, GIFT_db3$Code_function)],
    TRUE ~ trait
  )) %>%
  mutate(trait=factor(trait,levels=unique(GIFT_db3$Function))) %>%
  ggplot(aes(x=value, y=time_point, group=time_point, fill=type, color=type)) +
  geom_boxplot() +
  scale_color_manual(name="type",
                     breaks=c("Control","Hot_control", "Treatment"),
                     labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_fill_manual(name="type",
                    breaks=c("Control","Hot_control", "Treatment"),
                    labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
                    values=c("#4477AA50","#d57d2c50","#76b18350")) +
  facet_grid(trait ~ type, space="free", scales="free") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        strip.text.y = element_text(angle = 0)) + 
  labs(y="Traits",x="Metabolic capacity index")

GIFTs_elements_community_merged<-GIFTs_elements_community %>%
    as.data.frame() %>%
    rownames_to_column(var="sample") %>%
    filter(sample!="AD69") %>%
    pivot_longer(!sample,names_to="trait",values_to="gift") %>%
    left_join(sample_metadata, by = join_by(sample == Tube_code))%>%
    filter(time_point=="1_Acclimation"|time_point == "6_Post-FMT2")%>%
    mutate(functionid = substr(trait, 1, 3)) %>%
    mutate(trait = case_when(
      trait %in% GIFT_db3$Code_element ~ GIFT_db3$Element[match(trait, GIFT_db3$Code_element)],
      TRUE ~ trait
    )) %>%
    mutate(functionid = case_when(
      functionid %in% GIFT_db3$Code_function ~ GIFT_db3$Function[match(functionid, GIFT_db3$Code_function)],
      TRUE ~ functionid
    )) %>%
    mutate(trait=factor(trait,levels=unique(GIFT_db3$Element))) %>%
    mutate(functionid=factor(functionid,levels=unique(GIFT_db3$Function)))

# Create an interaction variable for time_point and sample
GIFTs_elements_community_merged$interaction_var <- interaction(GIFTs_elements_community_merged$sample, GIFTs_elements_community_merged$time_point)
  
ggplot(GIFTs_elements_community_merged,aes(x=interaction_var,y=trait,fill=gift)) +
        geom_tile(colour="white", linewidth=0.2)+
        scale_fill_gradientn(colours=rev(c("#d53e4f", "#f46d43", "#fdae61", "#fee08b", "#e6f598", "#abdda4", "#ddf1da")))+
        facet_grid(functionid ~ type, scales="free",space="free") +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=5),
              strip.text.y = element_text(angle = 0)) + 
        labs(y="Traits",x="Time_point",fill="GIFT")+
  scale_x_discrete(labels = function(x) gsub(".*\\.", "", x))

8.7 Wilcoxon comparison

8.7.1 Community elements differences: in CC acclimation vs post2

difference_table_CC %>% 
  filter(group_color=="Acclimation")
  Elements Acclimation      Post2                      Function     Element Difference group_color
1    T0503   0.1830714 0.01253516 Vitamin transport_Niacin (B3) Niacin (B3)  0.1705362 Acclimation

8.7.2 Community elements differences: in CI acclimation vs post2

difference_table_CI %>% 
  filter(group_color=="Acclimation")

8.7.3 Community elements differences: in WC acclimation vs post2

difference_table_WC %>% 
  filter(group_color=="Acclimation")

8.7.4 Comparison of both population in wild samples

difference_table_WC %>% 
  filter(group_color=="Acclimation")
  Elements Acclimation      Post2                              Function          Element Difference group_color
1    T0211  0.07898304 0.03960845           Amino acid transport_Serine           Serine 0.03937459 Acclimation
2    T0705  0.03026609 0.01792140      Nucleic acid transport_Allantoin        Allantoin 0.01234469 Acclimation
3    T0401  0.08571001 0.03355773    Organic anion transport_Spermidine       Spermidine 0.05215228 Acclimation
4    T0406  0.06660075 0.02039130     Organic anion transport_Succinate        Succinate 0.04620945 Acclimation
5    D0301  0.03445823 0.02175758             Sugar degradation_Lactose          Lactose 0.01270065 Acclimation
6    B0711  0.31689890 0.23431920 Vitamin biosynthesis_Menaquinone (K2) Menaquinone (K2) 0.08257970 Acclimation
7    T0503  0.05906234 0.02697581         Vitamin transport_Niacin (B3)      Niacin (B3) 0.03208653 Acclimation

####Butiryc acid biosynthesis

8.7.5 Comparison of both population in acclimation samples

8.7.6 Comparison of CC and CI in post2 samples

8.7.7 Comparison of CI and WC in post2 samples